function triangulated_surface

end

function builddtm, ground, gridsize
	triangulate, ground[0,*], ground[1,*], triangles
	;qhull, ground, triangles, /delaunay
	dtm = trigrid(ground[0,*], ground[1,*], ground[2,*], triangles, NX=gridsize, NY=gridsize)
	return, dtm
end

function subtract_dtm, nong, dtm, dtmbounds
	gridsize = (size(dtm, /dimensions))[0] ;assuming square dtm/chm

	dtmxr = dtmbounds[1] - dtmbounds[0]
	dtmyr = dtmbounds[3] - dtmbounds[2]

	;scale all nonground points so they relevant height can be extracted from the dtm
	sclnong = [(nong[0,*] - dtmbounds[0]) / (dtmbounds[1] - dtmbounds[0]), (nong[1,*] - dtmbounds[2]) / $
			(dtmbounds[3] - dtmbounds[2]), (nong[2,*] - dtmbounds[4]) / (dtmbounds[5] - dtmbounds[4])]

	;Derive absolute height for each point
	absnong = nong
	absnong[2,*] = nong[2,*] - dtm[sclnong[0,*] * (gridsize - 1), sclnong[1,*] * (gridsize - 1)]
	return, absnong
end

;dtmbounds is 6 element array
;0,1 X min X max
;2,3 Y min Y max etc
function buildchm, nong, dtm, dtmbounds, getabsheight=getabsheight
	gridsize = (size(dtm, /dimensions))[0] ;assuming square dtm/chm

	dtmxr = dtmbounds[1] - dtmbounds[0]
	dtmyr = dtmbounds[3] - dtmbounds[2]

	;scale all nonground points so they relevant height can be extracted from the dtm
	sclnong = [(nong[0,*] - dtmbounds[0]) / (dtmbounds[1] - dtmbounds[0]), (nong[1,*] - dtmbounds[2]) / $
			(dtmbounds[3] - dtmbounds[2]), (nong[2,*] - dtmbounds[4]) / (dtmbounds[5] - dtmbounds[4])]

	;Derive absolute height for each point
	absnong = nong
	absnong[2,*] = nong[2,*] - dtm[sclnong[0,*] * (gridsize - 1), sclnong[1,*] * (gridsize - 1)]
	if (arg_present(getabsheight)) then getabsheight = absnong

	;dimensions
	mind = [min(absnong[0,*]), min(absnong[1,*]), min(absnong[2,*])]
	maxd = [max(absnong[0,*]), max(absnong[1,*]), max(absnong[2,*])]
	nongxr = maxd[0] - mind[0]
	nongyr = maxd[1] - mind[1]

	;allocate chmbins - not the final CHM, just for finding points, will be trigridded
	chmbins = lonarr(gridsize * (nongxr / dtmxr), gridsize * (nongyr / dtmyr))
	cbsize = size(chmbins, /dimensions)
	sclabsnong = [(nong[0,*] - mind[0]) / (maxd[0] - mind[0]), (nong[1,*] - mind[1]) / $
			(maxd[1] - mind[1])]

	;now examine each point that falls in a particular
	chmbins[*] = 0
	chmbins[sclabsnong[0,*] * (cbsize[0] - 1), sclabsnong[1,*] * (cbsize[1] - 1)] = chmbins[sclabsnong[0,*] * (cbsize[0] - 1), sclabsnong[1,*] * (cbsize[1] - 1)] > absnong[2,*]

	arrin = array_indices(chmbins, where(chmbins gt 0))
	triangulate, arrin[0,*], arrin[1,*], arrtri
	chm = trigrid(arrin[0,*], arrin[1,*], chmbins[arrin[0,*], arrin[1,*]], arrtri, NX=cbsize[0], NY=cbsize[1])
	return, chm
end